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Abstract. We present a new radiative transfer code for axi-symmetric stellar atmospheres and compare test 
results against ID and 2D models with and without velocity fields. The code uses the short characteristic method 
with modifications to handle axi-symmetric and non-monotonic 3D wind velocities, and allows for distributed 
calculations. The formal solution along a characteristic is evaluated with a resolution that is proportional to 
the velocity gradient along the characteristic. This allows us to accurately map the variation of the opacities 
and emissivities as a function of frequency and spatial coordinates, but avoids unnecessary work in low velocity 
regions. We represent a characteristic with an impact-parameter vector p (a vector that is normal to the plane 
containing the characteristic and the origin) rather than the traditional unit vector in the direction of the ray. 
The code calculates the incoming intensities for the characteristics by a single latitudinal interpolation without 
any further interpolation in the radiation angles. Using this representation also provides a venue for distributed 
calculations since the radiative transfer can be done independently for each p. 
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1. Introduction 

Massive stars and their winds play an important role in 
shaping the dynamical structure and energy budget of 
galaxies. For example, they enrich the ISM with nuclear 
processed material and deposit large amounts of mechan- 
ical energy into their surroundings. Despite decades of 
research and considerable advancements in our under- 
standing of stellar envelopes, there is still much to learn. 
Because of the complexities of these systems, and the 
increasing emphasis on the details, it has become very dif- 
ficult to proceed without complex numerical simulations. 
It is not surprising, therefore, that the history of stellar 
studies reflects not only our advancing knowledge but also 
our increasing computational capabilities. Initially, simple 
plane-parallel LTE models were utilized in numerical 
simulations (see e.g.. lKurucJl99ll and references therein) 
and these were adequate for stars with dense atmospheres 
and low mass-loss rates. These models were also the 
only simulations that were viable on the computing 
facilities of the time. Unfortunately, the above simpli- 
fications cannot be extended to most early- type stars. 
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lAuer fc Mihaia^ lll972L HqtI . for example, demonstrated 
that the assumption of LTE is invalid in 0-type stars and 
the statistical equilibrium equations need to be solved for 
the level populations. For massive stars with extensive 
mass- loss (e.g., Wolf-Rayet stars) geometrical effects are 
also important and plane-parallel models are no longer 
sufficient. As a minimum, therefore, one needs to use 
non-LTE spherical models to understand these objects. 
The system of statistical equilibrium equations, however, 
is highly non-linear in the level populations and finding 
a solution for fully line blanketed models is a formidable 
task. We have reached the necessary level in computing 
power only in the last few years to be able to routinely 
perfo rm such computations (see e.g., 'Hubenv & Lan^ 
1995t ilauschil dt ct al.__ 1996.;_iIiUicr fc Miller. .1998; 
Pauldrach eTallboOlUGrlifener Vt alJl2002() . 

Plane-parallel and spherical non-LTE modeling have 
found wide applicability in spectroscopic studies. Recent 
works by ^Martins ct al. (2002); Crowthcr ct al. {20Q^; 
iHillier et alTij2003lh : lHerrero et alJ |2002^ have revised the 
temperature scale for O stars, for example, and have given 
new insights into the structure of stellar winds. However, 
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tions and cannot be used to study many important stellar 
objects. 

It has been known for a long time that some circum- 
stellar envelopes are non-spherical — the most well-known 
examples are the envelopes of Be stars. The hydrogen 
emission and infrared excess of these stars are thought to 
be produced in a thin disk. The presence of these disks was 
inferred fro m both line m odeling and polarimetric stud- 
ies (;Poeckert fc Marlborough 19 78a.,h,), and has been con- 
firmed by interferometr ic observations ijStee et al.lll995l: 
lOui rrenbach ct al.' "1997^. Further more, recent M HD sim- 
ulations (, ^Ca ssinclli ct al. 2002 : u d-Doula &: Owoc ki 2003; 
lOwocki fc ud-Doulal l2"oo1|r argue for equatorial confine- 
ment by magnetic field for the origin of the disks. If a 
dynamically important magnetic field is present in Be en- 
velopes that in itself ensures at least a 2D nature of their 
wind. 

Other stellar problems for which ID models are inad- 
equate include rapidly rotating OB stars, binaries with 
colliding winds or accretion disks, pre-main sequence 
and young stars, stellar envelopes irradiated by exter- 
nal sources (e.g., massive stars near an AG N), and the 
collapsing core (Typc-II) supernovae (e.g., IWang et alJ 
[200 1; Kifonidis et al. 2003.). Advanced supernovae mod- 
els may even have cosmological applications since these 
luminous objects can be used as distance calibrators in 
the nearby universe (see iDessart k, Hillier 2005a b . and 
references therein). 

The case of rapid OB rotators is particularly impor- 
tant for this paper since we test our code on such a 
problem. These stars are subjects of intense research and 
the exact structure of the rotating envelope is not well 
established. The conservation of angular momentum in 
the wind may result in meridional flow toward the equa- 
tor which potentially leads to disk formation (see e.g., 
biorkman & CassineUi 1993). Conversely the latitudinal 
variation of the surface gravity will result in a variation of 
the radiative flux with latitude that can inhibit disk for- 
mation, and can cause a stron g polar wind ( Owo cki et_.a,L 
ll99(THMaeder fc Mevnetl200ni) . Either way, the underlying 
spherical symmetry of the outflow is broken and at least 
axi-symmetric models are needed for spectral analysis. 

Motivated by the need for 2D model atmospheres, and 
by the availability of fast computers and methods, we un- 
dertook a project to develop a tool for spectroscopic anal- 
ysis of axi-symmetric stellar envelopes. The solution of the 
statistical equilibrium equations for the level populations 
and temperatu re is discus sed in the first paper of this se- 
ries l|Georgiev et ahlfeoOfil Paper I). At present the main 
code, ASTAROTH, solves for the radiation field by a con- 
tinuum transfer r outin e that is based on the method of 
iBusche fc Hilliej 1)2000(1 and uses the Sobolev approxima- 
tion for line transfer. In this paper we present an alternate 
routine for ASTAROTH that can handle the line-transfer 
without the use of Sobolev approximation in models with 
continuous, but not necessarily monotonic, velocity fields. 
We treated this problem independently from the main 



solution methods. In iJ21we describe our goals and motiva- 
tions in finding the proper solution method, and we also 
give a brief discussion of the chosen approach. The C-I--I- 
code that was developed for the transfer is described in 
where we also present the test results and verification. 
Finally, we draw our conclusions in 21 

2. Description of the Solution Technique 

A non-LTE model of a stellar envelope is a complex non- 
linear problem. The level populations and the radiation 
field are strongly coupled. Thus, an iterative procedure is 
needed to achieve a consistent solution. To solve the statis- 
tical equilibrium equations for the level populations, one 
must determine the radiative transition rates for free-free, 
bound-free and bound-bound transitions. These require 
the knowledge of the radiation moments 

J{v,y)^^ jj{v,n,v)d^ (1) 

and 

li{v)^^( f I{r,n,iy)<i>i{i^) diydil . (2) 
47r Jn Jq 

The quantities /(r , n, i/) , r , and n are the specific intensity, 
the spatial position, and the direction in which the radi- 
ation is propagating, respectively. The function rep- 
resents the normalized line-profile for any given bound- 
bound transition and the integrations are over all solid 
angles and frequencies. 

Only J and Ji are needed to solve the statistical equi- 
librium equations, but they have to be updated every it- 
eration cycle. This introduces stringent requirements on 
numerical efficiency and speed, but also allows for simpli- 
fications. The Radiative Transfer (RT) code does not have 
to produce the observed spectrum, for example, since it is 
irrelevant for the transition rates. Nor do the specific in- 
tensities at each depth need to be stored. On the other 
hand, the run time characteristics of the code are critical 
for its application in an iterative procedure. Therefore, our 
RT code is optimized to calculate J, Ji , and the "approxi- 
mate lambda operator" (A* , see i)2.3|l as efficiently as pos- 
sible. Crude spectra in the observer's frame are calculated 
only if requested, and only for monitoring the behavior of 
the code. 

At a minimum, a realistic non-LTE and line-blanketed 
model atmosphere requires the inclusion of most H, He, 
C, N, O, and a large fraction of Fe transitions in the cal- 
culation. The running time and memory requirements of 
such a model can be several orders of magnitude larger in 
2D than those of its spherical or plane-parallel counter- 
part. The dramatic increase in computational effort arises 
from both the extra spatial dimension, and from the ex- 
tra variable needed to describe the angular variation of 
the radiation field. In spherical models, for example, the 
radiation field is symmetric around the radial direction 
— a symmetry which is lost in 2D. We believe that re- 
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non-monotonic flow velocities, will inevitably require the 
simultaneous use of multiple processors. Therefore, we de- 
veloped ASTAROTH and this RT code to be suitable for 
distributed calculations by ensuring that their sub-tasks 
are as independent from each other as possible. 

2.1. The Solution of the Radiative Transfer 

Our choice to calculate moments J and J/ is to solve the 
radiative transfer equation for static and non-relativistic 
media 

nV/(r,n, t/) = -X(r,n, I/) [/(r,n,j/) - 5'(r,n,j/)] , (3) 

and then evaluate the integrals in Eqs.n]and[5] The quan- 
tities X 3.nd S in Eq. |3| are the opacity and source func- 
tion, respectively. A major simplification in this approach 
is that a formal solution exists for Eq.|21 At any s position 
along a given ray (or characteristic), the optical depth and 
the specific intensity are 

ru= I Xds' (4) 
Jo 

and 

I{t,) = Ibc e-"" + r ^(^') ' (5) 

Jo 

respectively (from now on, we stop indicating functional 
dependence of quantities on r, n, and v). Therefore, the 
intensity can be calculated by specifying Ibc at the up- 
stream end (s — 0) of the ray and by evaluating two inte- 
grals (assuming that S and x are known) . We sample the 
radiation field by a number of rays for every spatial point. 
If the number and the orientation of the rays are chosen 
properly, then the angular variation of / is sufficiently re- 
produced and accurate J and J; can be calculated. 

There are alternatives to this simple approach; each 
has its own merits and drawbacks. For example, from 
Eq. 1^1 one can derive differential equations for the mo- 
ments of the radiation field and solve for them directly. 
This approach has been successfully used in ID codes, like 
CMFGEN l|Hillier fc MiUeJ IlQQSV' and in calculations for 
2D continuum/grey problems (Busche 2001). A distinct 
advantage of the method is that electron scattering (ES) 
is explicitly included in the equations, and consequently 
no ES iteration is needed. However, to achieve a closed sys- 
tem of moment equations a closure relationship between 
the various moments is required. This relationship is gen- 
erally derived from the formal solution which requires at 
least a fast and rudimentary evaluation of Eqs. 0] and [3 
Furthermore, the 2D moment equations are quite compli- 
cated and it is not easy to formulate the proper bound- 
ary conditions in the presence of non-monotonic velocity 
fields. For our purposes we needed a simple approach that 
is flexible enough to implement in distributed calculations. 

An increasingly popular method to solve the RT is 
using Monte-Carlo simulations. In this method, a large 
number of photon packets are followed through the enve- 
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Fig. 1. A sub-section of a typical spatial grid used in our 
RT code. The boundary and internal points are indicated 
by grey and black dots, respectively. The solid arrow rep- 
resents a SC belonging to point i-f 2 and pointing in the 
direction of the radiation. Note, that the characteristic is 
terminated at the closest cell boundary (between nodes 2 
and 3), and is not followed all the way to the boundary 
of the domain (grey points). The numbering at the nodes 
indicates the order in which the intensity in this direc- 
tion is evaluated. The small empty circles on the SC are 
the integration points (see ti2.1|) and the dashed arrows 
show which grid points are used for interpolating x and S 
(straight arrows), or Ibc (curved arrows). 

by using this photon ensemble (see e.g., 
2003). While the Monte-Carlo simulations are flexible and 
suitable for parallel computing, they can also have unde- 
sirable run-time characteristics. It is also unclear how line 
overlaps in the presence of a non-monotonic velocity field 
can be treated by Monte-Carlo techniques without the use 
of Sobolev approximation. 

After considering our needs and options, we decided to 
use the straightforward approach, solving Eq. |2|and eval- 
uating Eqs. Hand 121 This approach provides a reasonable 
compromise of accuracy, numerical efficiency, and flexibil- 
ity. Our code will also increase the pool of available RT 
programs in stellar studies. Each solution technique has its 
specific strength (e.g., our method is fast enough for an 
iterative procedure) and weaknesses; therefore, future re- 
searchers will have more options to choose the best method 
for their needs. Having a selection of RT codes that are 
based on different solution methods will also allow for ap- 
propriate cross-checking of newly developed programs. 

The most accurate solutions for Eqs. 01 and El are 
achieved when the integrals are evaluated all the way 
to the boundary of the modeling domai n along each 
ray (Long Characteristi c (LC) method, Ijonesi Il973l: 
I.Tones fc Skumani'ci]ll973l) . To increase efficiency, we de- 
cided to use the so-called "Short-Characteristic" (SC) 
method, first explored by iMihalas et al.. (jl978|] and 
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method, the characteristics are terminated at the next up- 
stream radial shcU (normaUy, they would be terminated 
at any cell boundary) where Ibc is calculated by an inter- 
polation between the specific intensities of the nearest lat- 
itudinal grid points (see Fig.0. We calculate the specific 
intensity in a given direction for all grid points starting 
with those at the upstream end of the domain (where / is 
set to the appropriate boundary condition) and proceed 
with the calculation downstream (see Fig. for details). 
This evaluation scheme ensures that all intensity values 
are calculated by the time they are needed for the interpo- 
lation of Ibc- With this simple trick, the specific intensity 
is calculated very efficiently but for the cost of introducing 
coupling between the directional sampling of the intensity 
at the grid points. We will discuss the implications of this 
coupling in ii2.2l 

On every SC, we evaluate the integrals of Eqs.0]and 
|S1 for every co-moving frequency of the down-stream end 
point {i+2 in FigP) by 

r = At, At, = ^l±l±2^(.,.+i - s,) (6) 
and 

Jo 

N-l 

j=i J 

where iV-1 is the number of integration steps. Eqs.Eland 
[7| can be easily derived from Eqs. ^ and [S] by assuming 
that in each interval x ^^nd S are linear in s and r, respec- 
tively. To ensure that the spatial and frequency variations 
of the opacity and source function are mapped properly, 
we divide the SC into small Sj+i — Sj intervals by plac- 
ing enough "integration" points on the characteristic. The 
number of these points (N) depends on the ratio of the 
"maximum line of sight velocity difference" along the SC 
and an adjustable "maximum allowed velocity difference" . 
By choosing this free parameter properly we ensure ade- 
quate frequency mapping but avoid unnecessary calcula- 
tions in low velocity regions. Further, we can trade ac- 
curacy for speed at the early stages of the iteration and 
later "slow down" for accuracy. We allowed for 20 km s^^ 
velocity differences along any SC in the calculations that 
we present here. Even though this is larger than the av- 
erage frequency resolution of our opacity and emissivity 
data ('^10 km s~^), it was still adequate. Trial runs with 
2 km and 20 km s~^ "maximum allowed velocity dif- 
ference" for the ID model with realistic wind velocities 
(see iiH.2|l produced nearly identical results. 

The line of sight velocities, Xjj and Sj are calculated 
at the integration points by bi-linear interpolations using 
the four closest spatial grid points (see Appendix and 
Fig. Pi- We would like to emphasize, that the interpolated 



in which the integration is performed. This difference must 
be taken into account in Eqs.EHJjby applying the proper 
Doppler shifts at each integration point (see AppendixEl ■ 
With the exception of the intensity, all quantities are 
interpolated assuming that they vary linearly between 
nodes. Extensive testing of our code revealed that at least 
a third-order interpolation is necessary to calculate Ibc 
sufficiently accurately (see Appendix IA.2|I . For all other 
quantities first-order approximation is adequate in most 
cases but not in all. Since we wished to keep the first-order 
approximations if possible (it is the least time consuming 
and is numerically well behaved), a simple multi-grid ap- 
proach was introduced to improve accuracy. Unlike the 
intensity calculation, the interpolation of x ^md S do not 
have to be performed on the main grid; therefore, a dense 
spatial grid for opacities and source functi ons can be cre - 
ated, using monotonic cubic interpolation ()Steffenlh99(i|) . 
before the start of the calculation. Then, we use this dense 
grid to perform the bi-linear interpolations to the integra- 
tion points but perform the RT calculation only for spa- 
tial points on the main grid. Before the next iteration, the 
opacities and source terms on the dense grid are updated. 
To ensure a straightforward A* calculation we require the 
main grid to be a sub-grid of the dense grid. Further, the 
use of the dense grid is optional and only required if more 
accurate approximations of x and S are desired. With this 
rudimentary multi-grid technique, we improved the accu- 
racy of our calculations for essentially no cost in running 
time (~5-10% increase). However, there was a substan- 
tial increase in memory requirement. To avoid depleting 
the available memory, the RT is usually performed in fre- 
quency batches that can be tailored to fit into the available 
memory. This technique not only decreases the memory 
requirements, but also provides an excellent opportunity 
for parallelization. 

2.2. Our Coordinate System and Representation of 
Directions 

Most 2D problems that we are going to treat are "near- 
spherical" with a moderate departure from a general 
spherical symmetry. The radiation field is usually domi- 
nated by a central source in these cases, and it is practical 
to treat them in a spherical coordinate system. Therefore, 
we decided to use r, (3, and e (see Figure El for definition) 
for reference in our code. 

In spherical symmetry, the most natural way to map 
the directional variation of the intensity is using the "so- 
called" radiation coordinates, 9 and (j>, that are defined 

by 

cos{9) = n-r (8) 
and 

sin{9) ■ sin{(3) ■ cos{4>) = [n x r] • [r x z] . (9) 
The unit vectors n, r, and z are pointing in the direction 
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Fig. 2. The definition of our fundamental coordinate sys- 
tem. The unit vector n describes a characteristic (long 
thin line) pointing in the direction of the radiation and 
r, /3, and e are the traditional polar coordinates of a spa- 
tial point. Note that it is assumed here and in the rest of 
the paper that z axis is the axis of symmetry. We use the 
impact-parameter vector p (which is perpendicular to the 
plane containing the characteristic and the origin), instead 
of n, to represent a particular characteristic (see ii2.2l for 
explanations). 




Fig. 3. Diagram illustrating the connection between the 
radiation angle and the inclination angle. The "plane of 
the radiation" includes the characteristic and the origin. 
Angles i and /3 are the angular distances between the z- 
axis and the directions of the p and r vectors, respectively. 
Eq. El can be derived by a spherical sine law using the 
boldface spherical triangle. 



side of the z axis, respectively (see Fig.lJJ. A proper choice 
of Q angle grid can be very useful in treating inherent 
discontinuities around the limb of the central star and 
the symmetries due to the forward-peaking nature of the 



As mentioned in t l2.1l a serious drawback of the SC 
method is the interdependency of the specific intensities at 
different grid points. Beside introducing systematic errors 
by the successive intensity interpolations, the SC method 
also couples the directional sampling of the radiation field 
on the grid. Our choice of directions at a grid point not 
only has to suit the needs of the particular point but also 
has to be able to provide suitable starting values (Ibc) for 
other points. Unfortunately, 6 and (f) vary along a charac- 
teristic so it is not possible to use a uniform 9 and (f> grid 
for all grid points without intensity interpolations in the 
radiation coordinates. The later option is not desirable 
for multidimensional RT. First, it requires a large amount 
of memory to store all intensities for the interpolation. 
Second, it makes the parallelization of the code difficult. 

To find a proper directional sampling method one 
needs to look for quantities that are conserved along a 
characteristic, like 

p = r X n , (10) 

which we call the "impact-parameter vector" (see Fig.|21). 
This vector describes all essential features of a character- 
istic and can be considered as an analog of the orbital mo- 
mentum vector in two body problems. Its absolute value 
p— IpI is the traditional impact-parameter and its orienta- 
tion defines the "orbital plane" of the radiation (the plane 
that contains the characteristic and the origin). Following 
this analogy one can define an "inclination" angle for this 
plane by 



p ■ cosii) = p ■ z 



(11) 



In our code we set up a universal grid in impact- 
parameters (p) and in inclination angles (i) for directional 
sampling. As opposed to the 9 and (p angles, the inclina- 
tion angle and the impact-parameter do not vary along a 
ray; therefore, intensities in the proper directions will be 
available for the interpolation when the transfer is solved 
for a given i and p. Using an impact-parameter grid to 
avoid interpolation in 6 angle has already been inco rpo- 
rated into previous works fe.g'.. lBusche fc Hillierl200nl) . By 
introducing the inclination angle grid we simply exploited 
the full potential of this approach. 

It is useful to examine the relationship between the 
radiation angles and our directional coordinates. The con- 
version is via 



sin(9) — — 
r 



and 



sin{(f)) 



cos{i) 
sin{l3) 



(12) 



(13) 



at each grid point. Equation 1131 can be easily derived 
by spherical trigonometry as illustrated by Fig. |3 One 
can see from Eqs. El and El that there is a degener- 
acy between "incoming" - "outgoing" , as well as between 
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rays are defined hy ^ < (p < |7r.) The radiation coor- 
dinates (0, (j)) and (tt — 0, TT — 0) are represented by the 
same (p, i) pair. Fortunately, the "switch-over" can only 
occur at certain spatial positions. For example, the incom- 
ing rays become outgoing only at r= p, so this is just a 
simple book-keeping problem. Nevertheless, one should al- 
ways bear this degeneracy in mind when doing the actual 
programming implementation of our method. 

There remains one important question. How exactly 
do we choose the actual impact-paramet er and inclination 
angle grid? We adopted the approach of iBusche fc Hillieil 
l|200fl) who used the radial grid and a number of "core 
rays" (p < rcore) for the impact-parameters. The core rays 
are added only if a central source with a radius Vcore is 
present in the model. This will provide a radius dependent 
sampling since only p < r can be used for a given r radius. 
Also, the sampling is uneven and sparser around 9= ^ 
than around 9= or tt. Nevertheless, this grid was proven 
to be adequate for near spherical problems and also very 
convenient to use. For example, it ensures that p= r (the 
switch-over from "incoming" to "outgoing" ray) is always 
a grid point. Similarly, we based our inclination angle grid 
on the P grid, although, we have the option to define it 
independently. If needed, extra inclination angles can also 
be included around i= ^ to increase the (f) angle resolution 
at higher latitude. 

Fig. 01 illustrates a typical inclination angle grid and 
the (jj-angle sampling it provides. For illustration purposes 
we use a hypothetical P grid of ^tt (equator), -^n (72°), 
^TT (54°), ^77 (36°), ^TT (18°), and (pole). Then, one 
may choose these P values and their corresponding com- 
plementary angles (tt-/?) for the inclination angle grid. By 
our definition, angles « < f sample the < ^ < tt range, 
while « > f covers the rest of the </> space (see Fig. 
The behavior of the (/)-angle sampling created by this in- 
clination angle grid is very similar to that of the ^-angle 
sampling provided by the radial grid. One can easily see 
from Eq. El that for a given /3 any i < ^ — /? has no so- 
lution for (j). The equatorial regions (/? ~ therefore, 
are well sampled in angle while there is only one valid 
inclination angle at /3= {i= This is reasonable in 
axi-symmetrical models, as long as the polar direction is 
also the axis of symmetry (as we explicitly assume) . The 
(ji-angle sampling is also uneven. The regions around (/)— 
and TT (local meridian) are better resolved than those 
around (j)= ^ and (j)= ^n. In i)3.1H3.3l we will demonstrate 
that our sampling method not only eliminates the need for 
interpolations in 9 and (j) angles, but sufficiently recovers 
the directional variation of the radiation at every point 
and is adequate for RT calculations in axi-symmetric en- 
velopes. 

2.3. Approximate Lambda Iteration 

A seemingly natural choice for the iteration between the 
RT and level populations is the notorious "A-iteration" . 



cycle are used to calculate new J and J; which in turn 
are used to update the populations. Unfortunately, this 
simple procedure fails to converge for large optical depths. 
Convergence is ensured, howeve r, by using the Accelerated 
Lambda Iteration (ALI; see e.g.. lRvbicki fc Hummeill99ll: 
Hubenv.,1992.) which takes some of the inherent coupling 
into account implicitly. The relationship between J and 
the source function S can be summarized as 

J = A[S] , (14) 

where the A operator can be derived from Eqs. ^ and 
O Both A operator and S depend on the level popu- 
lations, however, we can "precondition" A (i.e., use the 
pop ulations from the previou s cycle to evaluate it, see 
e.g-l Rvbicki fc Hummei]ll99l|) and only take the coupling 
through S into account to accelerate the iteration. In 2D, 
A in its entirety is too complicated to construct and time 
consuming to invert, which is necessary to take the cou- 
pling into account. We can, however, split the A oper- 
ator into an "easy-to-invert" A* (Approximate Lambda 
Operator) and the remaining "difficult" part by 

J = A*[S] +(A-A*)[S'] . (15) 

Then, we can precondition the "difficult" part by using 
the old populations, and accelerate the iteration by in- 
verting A*. Note, that the full A operator never needs to 
be constructed, only A* since 

(A - A*) [S'-^] = .P-^ - A* [S'-^] (16) 

where J*~^ and S^~^ is the moment and source term from 
the previous iteration cycle. 

The actual form of A* is a matter of choice as long 
as it can be easily inverted. The most practical in 2D is 
separating out the local contribution (i.e., diagonal part 
of the A operator when written in a matrix form). This 
is easy to calculate and has reasonably good convergence 
characteristics. During the evaluation of moments J and Ji 
(see i i2.1(l . we also calculate the diagonal A* operator. This 
is a fairly straightforward book-keeping since we just have 
to add up the weights used for the local source function 
during the integration of Eq. |S1 

We used the A* operator to accelerate the ES iterations 
in our test calculations (see the following sections) . Apart 
from the initial "hiccups" of code development, the opera- 
tor always worked as expected and produced the publishe d 
convergence characteristics l|Rvbicki fc Hummed Il99l|) . 
The implementation of the ALO iteration into the solu- 
tion of th e statistical equilibrium equation is discussed in 
IPaoer 3 . 

3. Code Verification and Test Results 

We have developed a C++ code that implements the so- 
lution technique described in Sj21 As mentioned in ^2.1\ 
we used a modified version of the traditional SC method 
by terminating the characteristics at the closest spherical 
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Fig. 4. The (;6-anglc plane at different latitudes as viewed by an observer facing the central star/object. The unit 
vector pointing out of the page is toward the observer. Each figure is centered on the line of sight of the observer 
and the equator is toward the bottom of the page. The figure was created for inclination angles of 0°, 18°, 36°, 54°, 
72°, 90°, 108°, 126°, 144°, 162°, and 180° which are indicated near the head of the arrows. The radiation angle 4> is 
measured counter-clockwise from the direction toward the equator as indicated on the outer rim of the circles. Panels 
a and b are for f3= 36° and 54°, respectively. For clarity, we assumed that the impact-parameter (p) of the rays is 
equal to r; therefore, any direction that we sample lies in the (f> plane. The figure shows that the <f)-angle coverage is 
latitude dependent and unevenly spaced. Note, for example, the absence of i— 0°, 18°, 36° (and their complementary 
angles) for /?= 36°. 



cell boundaries in direction). This modification allows 
us to avoid intensity interpolations in the radial direc- 
tion which increases the accuracy when a strong central 
source dominates the radiation field. The transfer calcu- 
lation for an impact-parameter (p) and inclination angle 
(i) pair is performed on an axi-symmetric torus with an 
opening-angle of 2i and which is truncated at the inner 
radius of r= max(p, Vcore)- This torus contains all spa- 
tial regions that a ray described by p and i can reach. The 
calculation starts at the outermost radius and proceeds in- 
ward, shell by shell, until the truncation radius is reached; 
then, the outgoing radiation is calculated in a similar man- 
ner by proceeding outward. At the outer boundary we set 
the incoming intensity to zero while either a diffusion ap- 
proximation or a Schuster-type boundary condition can 
be used at the truncation radius if it is equal to Vcore- In 
its present form, the code assumes top-bottom symmetry, 
however, this approximation can easily be relaxed to ac- 
commodate general axi-symmetric models. The RT calcu- 
lation for each {p, i) pair is independent from any other. 
The only information they share are the hydrodynamic 
structure of the envelope, the opacities, and emissivities; 
all of which can be provided by ASTAROTH. 

There are at least two major venues to accommodate 
multi-processor calculations in the code. One way is to 
distribute the (p, i) pairs among the available processors. 
To optimize the calculation one needs to resolve a non- 



grid points involved in the RT is not the same for all (p, i) 
pairs, so the duration of these calculations is not uniform. 
For example, the transfer for p= and i= ^ involves all 
spatial grid points, while the one for p= and i= in- 
cludes only the points lying in the equator. To use the full 
capacity of all processors at all times, a proper distribu- 
tion mechanism needs to be developed that allows for the 
differences between processors and the differences between 
{p, i) pairs. 

We also have the option to distribute the work among 
the processors by distributing the frequencies for which 
the RT is calculated. In this case, the work-load scales lin- 
early with the number of frequencies, so the distribution 
is straightforward. However, the lack of sufficient memory 
may prevent the distribution of all opacities and emissiv- 
ities and the processors may have information only over 
their own frequency range. To take the effects of velocity 
field into account at the limiting frequencies, we introduce 
overlaps b(rtw(x;n the frequency regions. 

So far, we have performed multi-machine calculations 
where the (p, i) pairs or frequency ranges were distributed 
by hand. The results of the distributed calculations were 
identical to those performed on a single machine. Work is 
under way to fully implement distributed calculations by 
using MPI protocols. Since our goal is to run the entire 
stellar atmosphere code on multiple processors, we will 
discuss the details of parallelization in a subsequent paper 
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In the following we describe the results of some basic 
tests of our code. First, we calculate the radiation field 
in static 2D problems with and without ES. Then, we 
present our results for realistic spherical problems with 
substantial wind velocities. Finally, we introduce rotation 
in a spherical model and demonstrate the ability of our 
code to handle 2D velocity fields. 

3.1. Static 2D Models 

The basic characteristics of our code were tested by per- 
forming simple calculations, ID and 2D models without 
velocity field. We used the re sults of a LC program devel- 
oped by iHillied l)l994 Il996l) as a benchmark. This code 
was extensively tested and verified by reproducing one di- 
mensional models as well as analytical solutions available 
for op tically thin stellar envelopes fe.g.- lBrown fc McLeanI 
Il977|) . It was also tested against Monte-Carlo simulations 
of more complicated models. 

Our code reproduced the results of the LC program 
within a few percent for all spherical and axi-symmetric 
models. It was proven to be very stable and was able to 
handle extreme cases with large optical depths. The most 
stringent tests were the transfer calculations in purely 
scattering atmospheres. In such cases, the necessary iter- 
ations accumulate the systematic errors which highlights 
any weakness in the program. Several ID and 2D scatter- 
ing models were run with ES optical depths varying be- 
tween 1 and 100. Figures El and El compare our results to 
those of the LC code for a model with electron scattering 
opacity distribution of 



Xe 



= 10- 



(17) 



No other source of opacity and emissivity was present in 
the model. At the stellar surface we employed a Schuster- 
type boundary condition of Ibc = Ij while Ibc= was 
used at the outer boundary. The ES iteration was termi- 
nated when < 0.001% had been achieved. This model 
is an ideal test case since the ES optical depth is large 
enough to require a substantial number of iterations to 
converge, but the convergence is fast enough to allow for 
experimenting with different spatial resolutions. 

For the results we present in Fig El the LC code was 
run with 60 radial and 11 latitudinal grid points. The 4> 
radiation angle was sampled in 11 directions evenly dis- 
tributed between and tt. This code assumes top-bottom 
and left-right symmetry around the equator (/3= ^) and 
the local meridian [(f)— 0), respectively, so only half of the 
f3 and space had to be sampled. The radial grid, supple- 
mented by 14 core rays, was used to map the radiation 
angle dependence (see ^12.21 for description) . We used a 
slightly modified radial and latitudinal grid in our code. 
We added 3 extra radial points between the 2 innermost 
depths of the original grid, and 6 extra latitudinal points 
were placed between (3= and 0.15 tt. These modifications 
substantially improved the transfer calculation deep in the 




Fig. 5. The percentage difference between the J moments 
calculated by our (J^c) and by the LC program (J/c) as a 
function of the depth index (0 and 59 are the indices of the 
outer-most and inner-most radial grid points, respectively) 
for different latitudes. The ES opacity in this model is 
described by Eq. El The symbols *, o, x, □, and A 
indicate the differences for (3= 0, O.Itt, 0.27r, O.Stt, 0.47r, 
and § , respectively. Our code systematically overestimates 
J in the outer regions (0-40) which is mostly due to the 
second order accuracy of the radial interpolations. Errors 
from other sources (e.g., latitudinal resolution, (p angle 
sampling) are most important at high-latitudes {[3 0.1- 
0.2 tt) but still contribute less than ^1%. 



of the angle was identical to that of the LC code. We 
based our inclination angle grid on the /3 grid and added 
4 extra inclination angles around ^ to improve the cov- 
erage at high latitudes. This grid resulted in a latitude 
dependent (p angle sampling. At the pole, the radiation 
was sampled in only 2 directions while on the equator 60 
angles between and 27r were used. Note, that our code 
does not assume left-right symmetry! 

Figures El and El show that we were able to reproduce 
the results of the LC code within ~2% accuracy, and the 
total radial flux is conserved within 1% level. It is also 
obvious that our code needs higher spatial resolution to 
achieve the accuracy of the LC code. This is expected 
since the LC program uses higher order approximations 
and adds extra spatial points when needed to increase the 
overall accuracy. In fact, one should not call the LC code 
a pure Nr= 60, Ni}= 11 model. The auxiliary points in- 
creased the real resolution. It is not surprising, on the 
other hand, that our code runs substantially faster on the 
same machine. The difference was between a factor of 10 
and 2, depending on the number of iterations needed to 
converge. Unfortunately, we have not yet introduced so- 
phist icated ac celeration techniques, like the Ng accelera- 
tion (|N3[l92i, so our code is not the most efficient when 
a very large number of iterations is needed. 

The agreement between our code and the LC code 
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Fig. 6. The percentage loss/gain of the total radial flux 
{H — J^^ HrT^dVL) with respect to the total flux emanat- 
ing from the stellar surface (Hcore) as a function of depth 
index for the 2D model described by Eq. El The flux is 
conserved within ~1%. 



creased. Satisfactory agreement could be achieved, how- 
ever, by increasing the radial resolution. Our test prob- 
lems and most of the real problems that we will address 
later are near spherical with a modest latitudinal varia- 
tion. The intensity reflects the strong radial dependency 
and, therefore, the radial resolution controls the overall 
accuracy. Fig. |S1 reveals another feature of our method 
that affects the accuracy. Our result is sensitive to the 
high-latitude behavior of the intensity for a given incli- 
nation angle and impact-parameter. At the high-latitude 
regions, a given inclination angle samples directions that 
can be almost parallel with the equator. Slightly differ- 
ent directions that are almost parallel with the equator 
can sample very different radiation in some axi-symmetric 
models, such as models with thin disks. Aggravating this 
problem, our method also uses fewer directions to map the 
radiation field at these high latitudes, unless extra incli- 
nation angles around ^ are included. This explains why 
we had to use extra latitudes and inclination angles to 
produce the result for Figs. El and El We would like to em- 
phasize, however, that these problems are important only 
in extreme axi-symmetric models (e.g, very thin disks or 
strong polar jets). Many times, as it will be demonstrated 
in the next sections, reasonable accuracies can be achieved 
on ordinary and simple grids. 

During the static 2D tests, we also experimented with 
the multi-grid capability of our code and verified its scal- 
ing behavior. Tests with progressively increasing spatial 
resolution showed that our code has second order accu- 
racy. By doubling the number of radial grid points, for 
example, the errors decreased roughly 4-fold. We also per- 
formed the ES iterations in multiple steps and at pro- 
gressively increasing resolution. First, a coarse grid was 
created (e.g., half of the nominal resolution) for a crude 



Star AV 83 

Sp. Type 07 laf 

log g 3.25 

R 19.6 Ro 

Te// 34000 K 
M 2.5x10"^ Mq yr-i 

Voo 900 km s"^ 

13" 2 



" - Power for CAK velocity law JCastor et al.lll975^ . 

terms, a second iteration was performed on the nominal 
grid. This "double iteration" scheme was generally a factor 
of two faster than a single iteration on the nominal grid. 
This approach will be a promising venue for fast iterations 
in combination with other acceleration techniques. 

3.2. ID Test Cases with Realistic Wind Velocities 

After performing static 2D tests, we applied our code to 
realistic ID atmospheres. The primary purpose of these 
tests was to verify our handling of realistic velocity fields. 
We used a well k nown and tested ID s tellar atmosphere 
code, CMFGEN l)Hinier fc MilleJll998|) . for comparison. 
Observed spectra for a CMFGEN model are calculated 
independently by an a uxiliary routine, CMF_FLUX (see 
iBusche fc Hillieii I2OO5I for a description). We compared 
our simulated observed spectra to those of CMF_FLUX. 

We have an extensive library of CMFGEN models to 
choose a benchmark for our tests. We picked AV 83, a 
supergiant in the SMC (see Table [H which was involved 
in a recent study of O stars (Hillier et al. 2003}. Accurate 
rotationally broadened spectra wi th different viewing an- 
gles are also available for this star teusc he fc Hillier^l200,'^^ 
which we will use for comparison in US. 31 A detailed de- 
scr iption of the CMFG EN models for AV 83 can be found 
in iHillier et alj l|20n3l) . We chose their model v34_36C 
(see Table ^ to test our code. The radial grid with 52 
depth points was adopted from this model. The impact- 
parameter grid which samples the 9 radiation angle was 
defined by the radial grid augmented by 15 core rays (see 
tl2.2l for details). Our simulation was run as a real 2D case 
with two latitudinal angles {(3= and ^). We used 3 incli- 
nation angles which resulted in transfer calculations for 2 
and 4 angles in the polar and the equatorial directions, 
respectively. The RT calculations were performed on fre- 
quency regions centered around strategic lines, like Ha. A 
coarse grid {Nr= 26, Np— 2) and the nominal (iV^^ 52, 
Nf^= 2) grid was used for the ES iteration as in the cases 
of static models (see ii3.1|l . Note, that our model is not 
a fully consistent solution because we did not solve for 
the level populations. We simply used the output opac- 
ities and emissivities of the converged CMFGEN model 
and calculated the RT for it. 

Figure[3shows the normalized J moment as a function 
of wavelength for the CIV AA1548-1552 doublet and Ha 
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Fig. 7. The normalized J moment as a fmiction of wavelength around the C IV AA1548-1552 doublet (left column) 
and Ha (right column) at different locations in the envelope of AV 83 (the stellar model is described in Table The 
top row of figures shows J at Vr ^ foo, the middle at Vr ^ O.lwoo, while the bottom row displays J in the hydrostatic 
atmosphere {vr ^ 0). Note, that all spectra are in the co- moving frame. The solid (thin) and dash-dotted (thick) lines 
were calculated by CMFGEN and our code, respectively. Even though this model is spherical, our code treated it as 
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Fig. 8. The observed spectrum around the C IV AA1548-1552 doublet calculated by CMF_FLUX (thin sohd line) and 
by our code (thick dash-dotted line). Note, that these spectra are in the observers' frame. 



CMFGEN are in good agreement, except that we resolve 
narrow lines better. CMFGEN solves the moment equa- 
tion in the co-moving frame, starting at the largest fre- 
quency. This procedure introduces bleeding which broad- 
ens the sharp lines. Our results are not affected by this 
bleeding since we use the formal solution. 

Figure |H1 shows the observed spectrum for AV 83 in 
the observer's frame. As is the case with the J moment 
the agreement between our code and CMF_FLUX is ex- 
cellent. In this case CMF_FLUX does a better job but 
this is expected. Our code is primarily for providing J 
and J I for the solution of the rate equations while it pro- 
duces observed spectra only for testing. The main purpose 
of CMF_FLUX, on the other hand, is to produce highly 
accurate spectra in the observer's frame. 

We would like to emphasize that our code did not 
need higher spatial resolution to reproduce the results of 
CMFGEN/CMF_FLUX, as opposed to some cases pre- 
sented in ii3.1l The pure scattering models of ii3.1l were 
extreme examples and were hard to reproduce. The com- 
parison with CMFGEN proves that our code can handle 
realistic problems at a reasonable spatial resolution. 

3.3. Tests with a Rotating Envelope 

As a final test for our SC code we ran simulations of semi- 
realistic 2D atmospheres. These were created by introduc- 
ing rotation in otherwise ID models. AV 83 offers a good 
opportunity for such an experiment. It has a slowly accel- 
erating wind and low terminal velocity that enhances the 
importance of the rotational velocities. Also, its spectrum 
contains numerous photospheric and wind features which 
behave differently in the presenc e of rotation. 

Capitalizing on these features FBusche fc Hilheil 1 20051) 
used AV 83 to test their code for calculating observed 



study of the observable rotation effects. They utilized the 
LC method and a very dense directional sampling to cal- 
culate the observed spectra for an arbitrary viewing an- 
gle. This code serves the same purpose for ASTAROTH 
as CMF_FLUX does for CMFGEN; to calculate very ac- 
curate observed spectra for an already converged model. 
Since our code produces observed spectra only for test- 
ing purposes and error assessment, the comparison pro- 
vides on ly a consistency check between the two codes. 
Further, iBusche fc Hillierl l)2005() do not calculate radia- 
tion moments, so we could only examine whether our re- 
sults behave as expected with respect to the ID moments 
of CMFGEN. 

The rotation in the envelope of AV 83 was intro- 
duced by using the Wind Co i npressed Disk model (W CD , 
iBiorkman fc Cassinellll993(l . lBusche fc Hilheil l)2005(l ran 
several calculations to study the different aspects of rota- 
tion. We adopted only those t hat were used to study the 
Resonance Zone Effects fRZE. IPetrenz fc PuIsIITqQ^ . To 
isolate RZE-s, the latitudinal velocities were set to zero 
and the density was left unaffected by the rotation (i.e., it 
was spherical). The azimuthal velocity in such simplified 
WCD cases is described by 



V4> 



fcore ■ I r)\ 
Veq ■ • Sin (P) 



(18) 



rsee: lBiorkman fc Cassinellll993tlB"usche fc Hillieill2005|) . 

For the maximum rotational speed on the stellar surface 
(vp.n) we adopted 250 km s~^ following Busche & Hillicij 
l|2005|) . The radial velocity in the WCD theory is described 
by a CAK velocity law, so we used the same radial veloc- 
ities as in H3.2I 

We again adop ted the radial grid of model v36_34C 
l|Hillier et al.ll2003^ and used three (3 angles (0, f , and 
f ). In addition to these grids we had a dense radial and 
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of opacities, emissivities, and velocities; and a coarse grid 
{Nr^ 26, Nf3= 2) for the ES iteration. We used 14 in- 
clination angles, evenly spaced between and tt, which 
resulted in intensity calculations for 24 (p angles (between 
and 27r) at every point on the equator. As before, we per- 
formed our "double ES iteration scheme" (see S ii3.1l and 
I3.2|l with convergence criteria of < 0.001%. 

Figures ini and El show the behavior of the J moment 
around the CIV AA1548-1552 doublet and Ha, respec- 
tively, and also for the closest spherical and non-rotating 
CMFGEN model (thin/red line). It is obvious that sub- 
stantial deviation occurs only in the outer envelope and 
only for photospheric lines. Strong P-Cygni profiles, like 
those of the C IV AA1548-1552 doublet, are barely affected 
apart from a little smoothing around the blue absorption 
edge and at the maximum emission. The Ha emission, 
on the other hand, changes its strength substantially be- 
tween /3= and ^. This sensitivity casts doubts about 
the reliability of Ha as an accurate mass loss indicator for 
rotating stars with unknown viewing angle. A similar sen- 
sitivity to the rotation can also be seen on the iron lines 
around C IV AA1548-1552 doublet that are also formed at 
the wind base. 

Closer to the stellar surface the rotation effects on Ha 
diminish. At this depth, the behavior of narrow lines be- 
comes interesting. The iron lines around C IV AA1548- 
1552 doublet are broadened and skewed to the blue. This is 
the combined result of the large angular size of the stellar 
surface, limb darkening, and the broken forward-backward 
symmetry in the azimuthal direction. We will discuss this 
issue below in detail. At stellar surface {vr ~ 0) the optical 
depth is so large that any parcel of material sees only its 
immediate neighborhood which roughly moves with the 
same velocity. Consequently, no skewness, displacement 
or line-shape difference occurs between the profiles cal- 
culated for different latitudes (not shown in Figs. El and 

uni). 

Figure shows the detailed structure of the 
He I 4713.17 A profile in J moment. Since this line is not 
affected by blending (see e.g., the second panel of Fig. lT^ . 
its position, shape, and width should clearly reflect the ex- 
pected rotation effects and should highlight any inconsis- 
tencies in our model calculation. We present these profiles 
in velocity space and correct for the local radial veloc- 
ities. The bottom row of Fig. [Til shows He I 4713.17 A 
deep in the atmosphere {vr ^ and >> 1). The line is 
in weak emission centered around km as expected. 
The profiles are similar at all latitudes which reflects the 
fact that only radiation from the nearby co-moving regions 
contributes to J at this position. The line width reflects 
the local turbulent velocity and temperature. The top row 
of Fig. 111! shows the normalized J at ~ Voo- Here the 
line is in absorption and the profile widths show strong 
latitudinal dependence. We expect He I 4713.17 A to form 
in the photosphere, far from the radii where Vr ^ v^o 
(r ~ SOrcore)- In the co- moving frame of this position the 
central star covers only a small solid angle on the sky and 



ity, roughly equal to Woo- When we correct for the radial 
velocity of this position, we almost correctly account for 
the Doppler shift of each small section of the photosphere, 
hence, the profiles in Fig. ^2 should be and are centered 
on ~ km s^^. The polar view (solid line) shows the in- 
trinsic line profile (unaffected by rotation) while the equa- 
torial view (dash-dotted) broadened by ±250 km s^^ as 
it should. 

The profiles displayed in the middle panel of Fig. ITTI 
are more difficult to understand. They appear to be 
blueshifted and also skewed at (3= j and ^. At these in- 
termediate radii {vj. ~ O.Iwqo and r ~ 1.5rcore) the stellar 
surface covers a large portion of the sky and the Doppler 
shifts of photospheric regions vary substantially. The line 
profile in J is a superposition of the profiles emanating 
from different photospheric regions, and it is affected by 
the angular size of the photosphere and by the limb dark- 
ening. The line center should be rcdshiftcd by less than 
O.luoo velocity which explains the ^ —20 km s^^ blueshift 
in the middle panel of Fig llll fi.e.. we over compensated 
the Doppler shift). The blueward tilt of the profiles at 
[3 = J and ^ is caused by the forward-backward asymme- 
try around the rotational axis. The trailing and leading 
side of the photosphere contributes a broader and nar- 
rower profile, respectively, which causes the blueward tilt. 
We can conclude, therefore, that the gross characteristics 
of the He I 4713.17 A fine profiles in Fig.Ullreflect the ex- 
pected features at all depths and reveal no inconsistencies 
in our method. 

FigurelT^shows the observed spectra at different view- 
ing angles around selected transitions. We also show the 
calculations of CMF_FLUX for the corresponding spher- 
ical model. Not surprisingly, the observed spectra reveal 
the same characteristics as those of J moment at large 
radii. For our purposes, the most important feature of 
Figs. | 12 | is the remarka ble similarity to Figs. 4 and 5 of 
iBuschefc Hillieii l)2005|) . Despite the limited ability of our 
code to produce observed spectra. Fig. ^1 shows all the 
qualitative features of the synthetic observations. Most 
of the differences are due to our treatment of the ES. 
Our code does not redistribute the scattered radiation 
in frequency spac e which would produce smoother fea- 
tures like those of iBusche fc HillieJ l)2005|) . Note, that we 
run CMF_FLUX with coherent ES for proper compari- 
son; therefore, the spherical symmetric spectra also show 
sharper features. 

4. Summary 

We have implemented the short-characteristic method 
into a radiation transfer code that can handle axi- 
symmetric stellar models with realistic wind-flow veloc- 
ities. This routine will replace the continuum transfer 
plus Sobolev approximation approach that is currently 
used in our a xi-symme tric stellar atmosphere program 
(ASTAROTH, IPaoer f i. The new transfer code allows 
for non-monotonic wind-flow and, therefore, will enhance 
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Fig. 9. The normalized J moment as a function of wavelength around the C IV AA1548-1552 doublet at Vr ~ Woo 
(top) and at Vr ~ OTwoo (bottom). The wind velocity is described by a simplified version of the WCD model, for which 
the polar velocities and the density enhancements were turned off (see text for description). The azimuthal rotation 
was calculated by Eq.^jwith Veq= 250 km s~^. The thin (red) curve is the basic spherical symmetric model of AV 83 
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models for Be stars, OB rotators, binaries with colliding 
winds or accretion disks, pre-main sequence and young 
stars, and for collapsing core (Type-II) supernovae. 

The most important improvements of our approach are 
the sampling method that we introduced to map the direc- 
tional variation of the radiation, and the flexible approach 
to allow for non-monotonic velocity fields. We use a global 
grid in impact-parameters and in inclination angles (the 
angle between the equator and the plane containing the 
ray and the origin), and solve the transfer independently 
for every pair of these parameters. The code calculates 
the incoming intensities for the characteristics - a neces- 
sary feature of the short-characteristic method - by a sin- 
gle latitudinal interpolation. Our approach eliminates the 
need for further interpolations in the radiation angles. The 
effects of the wind-flow are taken into account by adapt- 
ing the resolution along the characteristics to the gradient 
of the flow velocity. This method ensures the proper fre- 
quency mapping of the opacities and emissivities where it 
is needed, but avoids performing unnecessary work else- 
where. Furthermore, it also provides flexibility in trading 
accuracy for speed. 

The code also allows for distributed calculations. The 
work-load can be shared between the processors by ei- 
ther distributing the impact-parameter - inclination angle 
pairs for which the transfer is calculated or by assigning 
different frequency ranges to the processors. 

We tested our code on static 1D/2D pure scattering 
problems. In all cases, it reproduced the reference re- 
sult with an error of a few percent. More complex tests 
on realistic stellar envelopes, with and without rotation, 
were also performed. Our code reproduced the results of a 
well-tested ID code fCMFGEN. HiUier fc Mi ller 1998). as 
well as the expected features in 2D rotating atmospheres. 
These tests demonstrated the feasibility and accuracy of 
our method. In a subsequent paper, we will describe the 
implementation of our code into ASTAROTH and present 
the results of fully self-consistent 2D simulations. 
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Appendix A: Interpolation Methods 



A.2. Interpolation of the Intensities 



A.l. Linear Interpolations 

We used bi-linear interpolations to calculate opacity, 
source function, and line of sight velocity at non-grid posi- 
tions in our modeling domain. The values were calculated 
by a weighted average of the corresponding quantities at 
the nearest grid points by 



4 

1=1 



1=1 

and 

4 

n - V = W; ■ (n ■ V;) 



(A.l) 



(A.2) 



(A.3) 



1=1 



The nearest grid points are described as 1= 1 {ri, (3i), 1= 2 
(^ii 02): 1= 3 (r2, /3i), 1= 4 (r2, P2)', which set the weights 
in Eas. rOlfO to 



Wl 



W3 



r2 



/3-/32 

ri - 7-2 Pi- P2 
r -ri 13-/32 
r2 -ri (3i- P2 



W2 



W4 = 



r -r2 (3 - (3i 



n -r2 (32- (3i 
r -ri f3 - Pi 
f 2 -ri P2- Pi 



(A.4) 
(A.5) 



Coordinates r (r2 > r > n) and P {P2 > P > Pi) are the 
coordinates of the general (non-grid) position. Note, that 
the frequency dependent quantities were interpolated in 
the co-moving frame! 

Since the integrals in Eqs.0]and|5lare evaluated in the 
co-moving frame of the down-stream end point of the char- 
acteristic, the opacities and source functions calculated by 
Eas. IA.l1 and IA.2l need to be properly Doppler corrected 
for the evaluation of the integrals. For a frequency 1^ in the 
co-moving frame of the down-stream end point the pro- 
cedure goes as follows: a; first we find the Doppler shifts 
Azj by Ea.rOlfor every integration point j on the char- 
acteristics (see t l2.1l for definitions) b; we find co-moving 
frequencies and Vk-i so that 



Linear interpolations of the intensities at the upstream 
end point of the SC does not provide acceptable accu- 
racy. This is because of the accumulation of errors from 
all previous intensity interpolations. Extensive testing of 
our method showed that the best result was achieved by 
using monotonic cubic interpolations fe.g. . .Steffen 199Q.) . 
We use this method to interpolate the intensities in P an- 
gle for fixed r and v. The monotonic cubic approximation 
provides the necessary 3rd order accuracy, yet avoids ar- 
tificial variations ("ringing") that can be amplified and 
propagated on our grid. Using monotonic interpolations 
actually dampens out such "ringings" and stabilizes our 
method. 

An unfortunate effect of requiring monotonicity, how- 
ever, is that we need to save the intensities for all frequen- 
cies and P angles on the previously treated shell (see ^^Ifor 
description of our code). Fortunately, this does not impede 
our efforts to accommodate multi-processor calculations 
because all {p,i) pairs can still be treated independently. 
Only, we require an additional memory area for NpNi, 
real number per (p, i) pair. 



Vk>v-{1- Azj) > Vk-i 



(A.6) 



at all integration points c; we find the opacity and source 
function for and Vk-i by Eas. IA.l1 and r^~2l d: we use 
linear interpolation in frequency space to get these param- 
eters at 1/ • (1 — Azj) for the integrations. 

This seemingly cumbersome procedure is actually a 
straightforward book-keeping that can be programmed 
very efficiently in the presence of monotonic velocity fields. 
Note, that we do not mean global monotonicity but a ve- 
locity field that is monotonic along the SC! One can assure 
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Fig. 10. Same as figure^ but for the spectra around Ha. 
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Fig. 11. Normalized J profiles of He I A4713.17 at Vr ~ Woo (top), O.luoo (middle), and (bottom). The solid, dashed, 
and dash-dotted spectra are for (3— (pole), and ^ (equator), respectively. The velocity scale is centered on the line 
and corrected for the above radial velocities. Our code reproduces the expected characteristics of the profiles within 
the uncertainties of our calculations (~ 20 km s~^). Note the skewed line profiles at intermediate radii (middle panel) 
which are the results of the broken forward-backward symmetry around the rotational axis (see ti3.3l for details). 
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Fig. 12. The observed spectra of AV 83 around the CIV AA1548-1552 doublet (top), the CIII/NIII/Hell emission 
complex between 4630-4700 A, Hell A5411, and Ha (bottom), respectively. See ti3.3l and Fig.|51for the description of 
the model parameters. The thick (red) curve is the spherical model calculated by CMF_FLUX, while the thin (blue), 
dashed (green), and dashed-dotted (purple) curves are our calculations for viewing angles 0, j, an d ^, respectively. 



